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Abstract - Using the momentum average approximation we show that the linear Holstein Hamil- 
tonian can be used to study single polarons only for weak electron-phonon couplings. For medium 
and strong linear coupling, quadratic and higher order electron-phonon coupling terms need to 
be included in the model. Doing so has drastic consequences for the properties of the polaron, 
which cannot be captured by a linear Holstein Hamiltonian with renormalized parameters. This 
shows that the linear Holstein Hamiltonian is inappropriate to model systems with strong electron- 
phonon coupling, at least for low carrier concentrations. 
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Introduction. — Coupling of carriers to phonons and 
the properties of the resulting quasiparticles, the polarons, 
are important for many materials, e.g. organic semicon- 
ductors [l], cuprates p], manganites [S], two-gap super- 
conductors like MgB2 [4], etc. In some cases the effective 
electron-phonon (el-ph) coupling A is known quite accu- 
rately. For others, like the cuprates, estimates range from 
very small (A ~ 0.3) to very large (A ~ 10) [sj. One pos- 
sible explanation for this is that, especially for stronger 
couplings where simple perturbational expressions are no 
longer valid, properly fitting the experimental data to the- 
oretical models can be quite involved [6]. 

Here we consider an even more serious issue: at strong 
el-ph coupling, the theoretical models themselves are in- 
consistent. In all widely-used models [7][8] one assumes 
at the outset that the displacements Xi of the atoms out 
of equilibrium are small enough to justify expanding the 
electron-lattice interactions to linear order in Xi. These 
linear models generically predict the formation of small 
polarons or bipolarons at strong coupling, with the car- 
rier(s) surrounded by a robust phonon cloud. As a result, 
lattice distortions {xi) can be considerable near the car- 
rier(s), invalidating the original assumption that higher 
order terms can be ignored in the Hamiltonian. 

In this Letter we investigate this issue in the single po- 
laron limits relevant for the study of weakly doped materi- 
als like very underdoped cuprates 9j and organic semicon- 
ductors [l], and for cold atoms/molecules simulators 10 . 



We study the ground-state (GS) of a single polaron in a 
generalized Holstein model including el-ph coupling up to 
quartic order in xi, to test the importance of the higher 
order terms. We find that for strong linear coupling even 
very small non-linear terms drastically change the proper- 
ties of the polaron. Moreover, we show that these effects 
go beyond a mere renormalization of the parameters of the 
linear Holstein model. As a result, attempts to find effec- 
tive parameters appropriate for a linear model by using its 
predictions to fit the properties of real systems are doomed 
to failure, as different values will be obtained from fitting 
different properties. This offers another possible explana- 
tion for the wide range of estimates of the el-ph coupling 
in some materials. More importantly, it means that we 
must seriously reconsider how to characterize such inter- 
actions when they are strong. Furthermore, this calls for 
similar investigations of the validity of these linear models 
at finite carrier concentrations, since it is reasonable to 
expect that they also fail in the strong coupling limit. 

To the best of our knowledge, we present here the first 
systematic, non-perturbative study of the importance of 
higher-order el-ph coupling terms on single polaron prop- 
erties. We note that in previous work going beyond linear 
models, purely quadratic (no linear term) but weak el-ph 
coupling was discussed for organic metals using pertur- 
bation theory 11 . In the context of high-Tc supercon- 



ductivity, the effect of anharmonic el-ph coupling on the 
phonons was studied in Rcf. 12 1, where a subset of the 
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diagrams from a quadratic electron-phonon coupling are 
included to obtain an effective coupling A. 

Formalism. — We use the momentum average (MA) 
approximation to carry out this study. MA was shown 
to be very accurate in describing GS polaron properties 
for the linear Holstein model, where it satisfies exactly 
multiple sum rules and becomes asymptotically exact in 
the limit of strong coupling |13|. It is straightforward to 
verify that all these considerations remain valid for the 
generalized Holstein model: 
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(1) 



where 'Hci = X^/e ^f'^k^k electron's kinetic en- 

ergy on a d-dimensional simple cubic lattice, with tk = 
"2^X^0=1 cos(fcQ,) for nearest-neighbor hopping (the car- 
rier's spin is irrelevant). — f^X^i^I^i describes a 



branch of Einstein phonons {% = 1). 



Finally, n^^_^^ 

9n^^c^^c^{b\ -\- bi)" describes hnear (n = 1) to quartic 
{n — 4) coupling terms. This general form is justified for 
models where the electron-lattice coupling modulates the 
on-site energy, of the type riiUixi), where = cjc^ is 
the density of electrons at a site, U{x) is the appropriate 
function describing the interaction between the additional 
electron and the atom it visits, and Xi b\ -\- bi is the 
displacement out of equilibrium of that atom. A Taylor 
expansion of U{x) up to quartic order leads to the gener- 
alized electron-phonon coupling of Eq. Q, where all the 
constants have been absorbed in the couplings <?„ which 
therefore have energy units. 

We define as the linear model the case where only gi ^ 
(i.e., the usual Holstein model); as the quadratic model the 
case where only gi 7^ 0, (72 7^ 0; and as the quartic model 
the case where all g„ ^ 0. The case with only 514 = is 
not considered because it is unstable since the potential 
felt by the ion where the electron is located has, in this 
case, the general form U{x) = ax^ + bx^ + ex. This has a 
minimum at — > 00, which is unphysical. 

The linear Holstein model is characterized by two di- 
mensionless parameters: the effective coupling strength 
A = gf/{2dtfl), where d is the dimension of the lattice, and 
the adiabaticity ratio n/(4(i<). As long as the latter is not 
very small, the former controls the phenomenology, with 
the crossover to small polaron physics occurring for A ^ 1 



15 . For ease of comparison, we continue to use these pa- 
rameters when characterizing the higher order models. For 
the quadratic model, the new energy scale 32 results in a 
third dimensionless parameter C = g2/gi- For the quartic 
model we could introduce two more parameters 173 / gi and 
(74/(71, but this makes the parameter space inconvenient 
to explore. To simplify it, we assume that gn = 5iC"~^- 
This choice is valid if the potential between the carrier and 
the ions is taken to be an unscreened Coulomb potential. 
Choosing a different potential U{x) will introduce some 



numerical prefactors in the definitions of gs, g^, which will 
modify slightly the results for the quartic Hamiltonian. 

We now describe in detail the MA solution for the 
quadratic model. The calculations for the quartic model 
are analogous but much more tedious. 

We want to find the single particle Green's function 
G(fc,w) = {0\c^G{cu)cl\0) where G(w) = [Lo-H + ir)]~^ is 
the resolvent for this Hamiltonian, with 77 — > a small 
positive number and |0) the vacuum state. From this 
we can extract all the polaron's GS properties 13 . We 
rewrite the quadratic Hamiltonian as "H = Hq - 
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Til, where 
while Hi = 



The equation of mo- 



92 cjc 

5i(6l+6.)+.92(for + 
tion (EOM) for the propagator is obtained recursively 
from Dyson's identity, G{uj) = Ga{^) + G'(a;)'HiG'o(w) 
where Go(a;) = [w — T^o + ^^y]"^ is the resolvent for 'Hq. 
Using it in G{k,u}) yields the EOM 



^ G(fc,^) = Go(fc,c^) 



EE 



g„F„(fc,w;i) 



(2) 



where Fn{k,uj;i) = (0|CfcG(a;)4(&J)"|0) . 

Applying Dyson's identity to generate EOM for the 
Fn propagators results in an infinite system of coupled 
equations which involves many other generalized propaga- 
tors. MA 13 , 14 circumvents this complication by mak- 
ing the approximation Go{i — j,u};n) « dijgo(uj]n) for any 
n > 1, where Go{i - j,Lj;n) = ^ (0|c,5f Go(t^)(6l)"4 10). 
This is justified because the polaron GS energy lies be- 
low the free particle spectrum, and for such energies 
the free-particle propagator decreases exponentially with 
|i — j|. Thus, MA keeps the largest contribution and ig- 
nores the exponentially smaller ones. This becomes ex- 
act in the strong-coupling limit t — > 0. The propaga- 
tor gQ{uj;n) — [l/gQ^u — nil — g2) ~ 2g2n\^^ is that of 
a carrier scattered by an on-site potential 2g2n, where 

ffo(w) = wJ2k'^/(^-^k + iv)- 

MA allows us to obtain a simplified hierarchy of EOM 
involving only the generalized Green's functions F„. For 
any n > 1, they read: 

Fn{k,uj;i) ^go{uj;n) ■ [n{n - l)g2Fn-2{k,(jJ]i) 
+ ngiFn-i{k,uj;i) + giFn+i{k,uj;i) + .92^+2 (fc,w; z)] . 

Since the arguments of all Fn propagators are the same, 
we suppress them in the following for simplicity. 
Following the technique introduced in Ref. 
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we re- 
duce this to a simple recursive relation for the vector 
W-„ = (F2„_i,F2„). The EOM for W„ are 7nW„ = 
anWn—i + l3nWn-\-i, whcrc the l3n and 7n are 2x2 
matrices whose coefficients are read off of the EOM, 
namely a„\ii = (2n - l)(2n - 2)g2go{uj;2n - 1), a„\i2 = 
(2n - l)gigo{ui; 2n - 1), a„\2i = and an|22 = 2n{2n - 
l)52ffo(w;2n), while 



fg29o{i^;2n - 1) 

V 5i5o(w;2n) g2go{uj\2n] 



(3) 
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(4) 



1.5 



~2ngi5o(w; 2n) 

This simple recursive relation for Wji has the solution 
Wn = An W„_i for any n > 1, where Ap are 2x2 matrices 
obtained from the infinite continued fraction 



[7n - ;5nAn+l]" 



In practice, we start with Ajf = for a sufficiently large 
cutoff N, chosen so that the results are insensitive to fur- 
ther increases in it (N ^ 100 is usually sufficient). 

,0 a22^ 

after using Eq. (|5| iV — 1 times. As a result, Fi ~ ai2i^o, 



, where ai2 and 022 are obtained 



We find Ai = 



F2 - 022^^0, where G(fc,cj) = J2^e'''--^ /VNFo{k,iu;i). 
Using these in Eq. ^ leads to a solution of the expected 
form G{k,bj) = [a; — — + irj] ^ , with the MA self- 
energy for the quadratic model: 



E(w) = 31012(0;) + 32022(0;). 



(6) 



The reason why the self-energy is local at this level of MA 
is the simplicity of this Hamiltonian, whose vertices are 
momentum independent; this issue is discussed at length 
for the linear Holstein model in Ref. [14 . 

The quartic model is solved analogously. The main dif- 
ference is that here the EOM for F„ involves 9 consecutive 
terms, from f„_4 to -^^+4. These can also be rewritten as 
simple recurrence relations 7nW„ = anWn—i + /3nWn-|-i, 
but now Qfn, /3n and "/„ are 4x4 matrices. Their expressions 
are too long to be listed here. 

Results and Discussion. — From the Green's func- 
tion we can find the polaron GS energy Eos , its quasipar- 
ticle weight Z (for local self-energies, the effective mass 
m* /m = 1/2'), the average number of phonons in the GS 
polaron cloud, iVp/j, etc. 

There are many options on how to gauge the importance 
of the non-linear terms. We begin with the most standard 
approach, namely by studying their effect on Egs- In the 
following, we use the rather ad hoc criterion that a higher 
order term must be included in the Hamiltonian if its in- 
clusion changes Egs by more than 10%. In our opinion, 
this is an upper estimate for an acceptable threshold; ob- 
viously, choosing a smaller value will make the non-linear 
terms even more relevant. In Fig. [l] we show the regions of 
validity for the linear and quadratic models, according to 
this rule, in the parameter space {(, A) for Q./t = 0.5. The 
solid line shows the boundary above which quadratic terms 
become important in ID (full and empty circles mark this 
boundary in 2D and 3D, respectively). The dashed line 
is the boundary above which quartic and higher order 
terms become important in ID, according to this crite- 
rion (full/empty squares mark the 2D/3D analogs). 

These results confirm the general expectation that while 
the use of linear models is safe for small A, the quadratic 
term becomes important with increasing A. What is sur- 
prising, however, is how fast the window of validity of the 
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Fig. 1: (Color online) Regions in the parameter space where 
using linear (n = 1), quadratic (n = 2) and higher order models 
(n > 4) is appropriate. Lines, full and empty symbols indicate 
the locations of the crossovers from n = 1 to n = 2 and n = 2 
to n = 4 in ID, 2D and 3D, respectively, for Q./t — 0.5. 



linear model closes down, even for this rather generous 
criterion. Indeed, Fig. [T] suggests that quadratic terms 
should be included for all intermediate and strong cou- 
plings A > 1. On the other hand, cubic and quartic terms 
are much less important for the physically acceptable val- 
ues C < Ij for ^-ny linear coupling A. 

The importance of the quadratic term at medium and 
large A is further validated and even amplified by its effects 
on other GS properties. Fig. [2] shows the quasiparticle 
weight Z — m/m* and the average phonon number Npfi 
vs. C for quadratic and quartic models in ID. Results in 
higher dimensions are qualitatively similar to these ID 
results for small A, and become quantitatively similar to 
them in the interesting regime of large A where all of them 
converge towards those of the atomic limit t = 0. 

First, we note that the C = intercepts trace the pre- 
dictions of the linear model: with increased coupling A, 
Z decreases while Nph increases as the polaron acquires 

From these intercepts, we 



a robust phonon cloud 13 15 



estimate that the linear model predicts the crossover to 
the small polaron regime to occur around A ^ 1.5 for this 
adiabaticity ratio and dimension. 

The quadratic model, whose predictions are indicated 
by lines, shows the dramatic importance of for strong 
linear coupling A > 1.5: here both Z and Nph vary by 
about an order of magnitude as increases from to 0.1. 
These changes are much more significant than those in 
Egs, since according to Fig. [T]the point A = 1.5, C = 0.1 
lies just above the threshold of 10% change in Egs, which 
is very modest compared to the order of magnitude vari- 
ation in Z and Nph- This shows that Fig. [T] should only 
be used as a rough guide for when the quadratic term 
becomes relevant, with the understanding that it may sig- 
nificantly underestimate its importance; indeed, for large 
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Fig. 2: GS quasiparticle weight (left panel) and GS average 
phonon number (right panel) vs. ^, in the quadratic (n — 2, 
lines) and quartic (n = 4, symbols) models, for various values 
of A and Q — 0.5t, in one dimension. 
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Fig. 3: (Color online) (left) Z^t, and (right) A^^^*' vs. (, for 
gi = -\/2 and = 0.5 (full lines). Dashed (red) lines show the 
mean-field estimates, while the dot-dashed (blue) lines show 
the results of fitting g/i} to exactly reproduce the other quan- 
tity. See text for more details. 



A this term cannot be ignored unless C ^ 0.1. 

The quadratic term completely changes the behavior 
of the polaron in the limit of medium and large A. For 
example, in the quadratic model at A = 1.5 and C ^ 0.1 
the polaron is light and with a small phonon cloud, in 
complete disagreement with the linear model prediction of 
a heavy small polaron at this A. For higher ^, Z and Np^ 
have a slight turnaround towards smaller/larger values, 
for reasons explained below, but are still consistent with 
a large polaron. Inclusion of cubic and quartic terms (the 
symbols in Fig. [2] show the results of the quartic model) 
further changes Z and Np/^, but these changes are much 
smaller for all C, of up to ^ 10% when compared to the 
quadratic model values, as opposed to order of magnitude 
changes between the quadratic and the linear models. This 
is in agreement with Fig. [ij which also suggests a very 
limited effect of the cubic and quartic terms for any C < 1- 
To understand the effects of the quadratic term at large 
A, we study it in the atomic limit t = {X = oo) where 
the carrier remains at one site and interacts only with the 
phonons of that site. Focusing on this site, its quadratic 
Hamiltonian "Wit^ = Qb^ + ELi5n(^^ + is well- 
studied in the field of quantum optics, where it describes 
so-called squeezed coherent states [Tt]. It is diagonalized 
by a displacement and a Bogoliubov transformation, with 
new bosonic operators 7^ = ub^ + vb + w, where u, v 



and w are such that Ti 



n 



(2) 



giy/n/nl^ and v = 



E^(^l\ We find 



= v/(^^ + 2g2+^^at)/(2nat), W = 
Sgn(g2)\/(f^ + 2^2 - f^at) /(2f7at). 



From these, we obtain Eq^ 
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and Za 



h exp [- 
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The latter result requires the expansion of the squeezed 
coherent states in the number state basis [Ts] . 

Figure [3] shows Z^t and A^^^*' C (thick lines), which 



agree well with the corresponding A = 2 results of Fig. [2] 
In particular, for ^ -> we find fl^t — + 2.giC + C(C^)i 
= ^ [1 - + C'(C^)] , explaining their linear in- 
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^ph 

crease/decreases for small (. 

The slight turnaround of the Z and iVph curves at larger 
values of ( is also observed in the atomic limit of the 
quadratic model. The reason is that the first term in N^^^ 
increases whereas the second term decreases with C. As 
discussed above, for small ( the second term dominates 
and the overall number of phonons decreases. For large 
however, the second term vanishes whereas the first term 
diverges as ^/g2 = VQh- Hence, as ( increases N^^ has a 
minimum, and then starts to increase with C. Basically, 
here the (72 (^^^ + b^) coupling dominates over the linear 
coupling gi{b^ + b) and changes the trend. 

This leads us to pose the question whether these ex- 
act results of the quadratic atomic model can be fit well 



by an effective linear model H 



(1) 



nb^b + g{b^ + b), for 



some appropriate choice of the effective parameters 51, g. 
One way to achieve this is with a mean-field ansatz k, 
2(6t)6t _ (6t)2^ ^ith the GS expectation value of 5^ 
The self-consistency condition {W) — —{gi+2g2{h^)) / {^ + 
2g2) leads to the mean-field estimates f^MF = 51-1- 2172, 
gjAF = 51 - 2gi52/(f^ + 4.g2). Thus, for small C = 52/51. 
JImf increases whereas ^mf decreases with increasing C 
so the effective coupling A — g^l{2dtVl) decreases with 
C^. This is consistent with the observed move away from 
the small polaron limit with increasing C,. Quantitatively, 
however, these mean-field results (dashed lines in Fig. [s]) 
are not very accurate for small C, and fail to capture even 
qualitatively the correct behavior when C ^ 1, since here 
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In fact, there is no choice for effective linear parame- 
ters g and fl that reproduces the results of the quadratic 
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model. This is because in the linear model, both Z and 
TVph are functions of gj^ only. Fig. [s] shows that if one 
chooses this ratio so that n'^^ = ^ph, then Z (dot-dashed 
line in the left panel) disagrees with Zat, and vice versa. 
Even more significant is the fact that even if one could find 
a way to choose ^ so that the overall agreement is sat- 
isfactory for all GS properties, the linear model's predic- 
tions for higher energy features would still be completely 
wrong. For example, it would predict the polaron-|-one- 
phonon continuum to occur at Eqs + i7 instead of the 
proper Eqs ~\- $7 threshold. Since in the atomic limit the 
predictions of the quadratic model cannot be reproduced 
with a renormalized linear model, we conclude that this 
must hold true at finite hopping t as well, at least for large 
A where the quadratic terms are important. 

So far we discussed moderate values of the adiabaticity 
ratio il/i = 0.5, as well as the anti-adiabatic (atomic) 
limit. MA predicts similar results in the adiabatic limit 
ri/t — > for large A, where it remains accurate, but is 
unsuitable to study small and moderate couplings [M] . We 
expect that here the quadratic coupling is essential even 
for small couphngs A — 0, because the term Igi ^ ■ 
ensures that phonons are gapped even though $1 = 0. 

The discussion so far was only for the case C > 0. The 
behavior of models with C < can be glimpsed at from 
the exact results in the atomic limit. For small negative 
C, the results listed above show that the average phonon 
number increases with |^| while the qp weight Zat de- 
creases fast, i.e. the polaron moves more strongly into the 
small polaron limit. This is in agreement with the MA 
predictions for the quadratic model (not shown). Here, 
however, we must limit ourselves to values |C| < ^^/(45i) 
so that r^at remains a real quantity (a similar threshold is 
found for the full quadratic model. Note that the value of 
this threshold decreases with increasing A). For values of 
ICI above this threshold the quadratic model becomes un- 
stable. This, of course, is unphysical. In reality, here one 
needs to include higher order (anharmonic) terms in the 
phonon Hamiltonian "Hph since they guarantee the stabil- 
ity of the lattice if the quadratic terms fail to do so. Such 
anharmonic terms may have little to no effect in the ab- 
sence of the carrier, but clearly become important in its 
presence, in this limit. They can be treated with the same 
MA formalism we used here. Their effects, as well as a full 
analysis of all possible signs of the non-linearities and the 
resulting polaron physics will be presented elsewhere. For 
our current purposes, it is obvious that in the case C < 0, 
higher order el-ph coupling terms also play a key role in 
determining the polaron properties unless A is very small, 
and therefore cannot be ignored. 

The results presented so far clearly demonstrate the im- 
portance of non-linear el-ph coupling terms if the linear 
coupling A is moderate or large, through their significant 
effects on the properties of a single Holstein polaron. A 
reasonable follow-up question is whether such dramatic 
effects arc limited to the single polaron limit or are ex- 
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Fig. 4: Estimate of the bipolaron phase diagram in ID for 
Q./t = 0.5 and for different values of C,, based on second or- 
der perturbation theory in t. In all four panels, the solid lines 
show the transition from So (on-site) stable bipolarons to Si 
(nearest-neighbor) stable bipolarons, while the dashed lines 
show the unbinding transition above which bound polarons 
are unstable. Note that panels (c) and (d) have a significantly 
rescaled y-axis. See text for more details. 



pected to extend to finite carrier concentrations. While 
the limit of large carrier concentrations remains to be in- 
vestigated in future work, here we present strong evidence 
that quadratic terms are likely to be equally important at 
small but finite carrier concentrations. 

Of course, for finite carrier concentrations one needs 
to supplement the Hamiltonian with a term describing 
carrier-carrier interactions. The simplest such term is an 
on-site Hubbard repulsion Hu = U J2i''^it''^ily ^^'^ gives 
rise to the Hubbard-Holstein Hamiltonian. The linear ver- 
sion of this Hamiltonian has been studied extensively by 

In particular, for low 



a variety of numerical methods 15 
carrier concentrations its phase diagram has been shown 
to consist of three regions: (i) for large A and small U, the 
deformation energy favors the formation of on-site bipo- 
larons, also known as the 5*0 bipolarons; (ii) increasing U 
eventually makes having two carriers at the same site too 
expensive, and the 5*0 bipolarons evolve into weakly-bound 
5*1 bipolarons, where the two carriers sit on neighboring 
sites. The binding is now provided by virtual hopping pro- 
cesses which allow each carrier to interact with the cloud 
of its neighbor. However, at small A and large U this 
binding mechanism is insufficient to stabilize the SI bipo- 
laron, and instead one finds (iii) a ground state consisting 
of unbound polarons. 

This phase diagram has been found numerically in ID 
and 2D 



19 



20 



for the linear Hubbard-Holstein model. 
Some results in 3D have also become available very re- 
cently 21 . In ID and 2D, the separation lines between 



the various phases are found to be close to those estimated 
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using second order perturbation theory in the hopping t, 
starting from the atomic Umit 19 20 . This is expected 



since for large hnear coupling A, the results always con- 
verge toward those predicted by the atomic limit. 

Since the quadratic Hamiltonian can be diagonalized 
exactly in the atomic limit, we use second order perturba- 
tion theory in t to estimate the location of the separation 
lines for various values of C > 0. The results are shown 
in Fig. [4] Panel (a) shows the rough phase diagram for 
C = 0, in agreement with the asymptotic estimates shown 
in Refs. [19^,^20^ (note that the definition of the effective 
coupling used in those works differs by various factors from 
our definition for A). Panels (b)-(d) show a very signifi- 
cant change with increasing Even the presence of an 
extremely small quadratic term ( — 0.01 moves the two 
lines to considerably lower U values, as shown in panel 
(b), while for C = 0.05 and 0.1, the bipolarons are stable 
only in a very narrow region with small values of U (note 
that the vertical axes are rescaled for panels (c) and (d)). 

The dramatic change with increasing C in the location of 
these asymptotic estimates for the various bipolaron tran- 
sitions/crossovers strongly suggests that non-linear el-ph 
coupling terms remain just as important in the limit of 
small carrier concentrations as they have been shown to be 
in the single polaron limit. In particular, these results sug- 
gest that the presence of non-linear el-ph coupling terms 
leads to a significant suppression of the phonon-mediated 
interaction between carriers, so that the addition of a small 
repulsion U suffices to break the bipolarons into unbound 
polarons (whose properties are also strongly affected by 
the non-linear terms, as already shown). 

The Holstein model is the simplest example of a g{q) 
model, i.e. a model where the electron-phonon interaction 
depends only on the momentum of the phonon. Physi- 
cally, such models appear when the coupling to the lattice 
manifests itself through a modulation of the on-site en- 
ergy of the carrier. The Frohlich model is another famous 
example of g{q) coupling. Models of this type are found 
to have qualitatively similar behavior, with small polarons 
forming when the effective coupling increases. These small 
polarons always have robust clouds, with significant dis- 
tortions of the lattice in their vicinity. We therefore expect 
that non-linear terms become important for all such mod- 
els at sufficiently large linear coupling. 

To summarize, we have shown that non-linear terms in 
the el-ph coupling must be included in a Holstein model 
if the linear coupling is large enough to predict small po- 
laron formation, and that doing so may very significantly 
change the results. We also argued that these changes 
cannot be accounted for by a linear Holstein model with 
renormalized parameters. These results show that we 
have to (re)consider carefully how we model interactions 
with phonons (more generally, with any bosons) in mate- 
rials where such interactions are expected to be strong, at 
least for low carrier concentrations and for models where 
this coupling modulates the on-site energy of the carriers. 
Whether this is also true in the metallic regime and/or for 



other types of models remains an open question. 
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